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ABSTRACT 


The GLAS model's surface fluxes of sensible and latent 
heat were found to exhibit strong 2>6t oscillations at the 
individual grid points as well as in the sonal and hemispheric 
averages. In addition, it was pointed out by Charney et al. 
{ 1911 ) that a basic weakness of the GLAS model has been its 
lower evaporation over oceans and higher evaporation over land 
in a typical monthly simulation. On examining the planetary 
boundary layer (PBL) parameterixation, it appeared that the 
calculation of surface temperature and the use of ^ ^qc 
constants in the eddy diffusivity calculation for the mixed 
layer were primarily responsible for these deficiencies. 

The GLAS model PBL parameterization has been changed 
to calculate the mixed-layer temperature gradient by solution 
of a guadratic equation for a stable PBL and by a curve-fit 
relation for an unstable PBL. The basic formulae used to 
determine the drag coefficient, its stability dependence, and 
the effect of moisture on the temperature gradient remain 
unchanged. The new PBL parameterization yields surface temper- 
ature and surface fluxes without any 2-6t oscillation. Also, 
the geographical distributions of the surface fluxes are im- 
proved. 

The parameterization presented here is incorporated 
into the new GLAS climate model. Seme results which compare the 
evaporation over land and ocean bttwoen old and now calculations 
arc appended. 
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INTRODUCTION 


In a series of climate simulations with the GlAS GCM, 
it was found by Charney et al« (1977) that the zonal) y averaged 
monthly mean evaporation was less than observed over the oceans 
and greater than observed over the land. This weakness in the 
model has persisted consistently in all summer and winter simu- 
lations. We also noticed a 2-6t oscillation in the evaporation. 
This oscillation was evident in the zonal and hemisphere averages 
as well as at individual grid points. At this point, a systematic 
examination of the boundary layer parameterization was undertaken. 
Before discussing the details of this problem, we provide the 
following background on the GLAS PBL parameterization. 

The present PBL parameterization was originally 
developed by Katayama^ (1972). Later, the Katayama parameteriza- 
tion was modified by Somerville et al. (1974), who introduced 
the formulation (Deardorff, 1967) for the eddy diffusivity of 
the mixed layer. For some reason, not discussed by Somerville 
et al. (1974), the original constants used by Deardorff were 
somewhat modified. In this parameterization, the PBL is assumed 
to be made up of two layers: the surface layer and the mixed 

layer. The surface variables are defined at the interface between 
these layers. These are surface temperature T^, the surface 
humidity qg, and surface wind components Ug and Vg. The 
surface layer is very shallow; its depth is about 10-50 m. For 
this reason, its heat and moisture capacities are negligible. 
Therefore, Tg is an equilibrium temperature determined by the 

^ As described by Arakawa (1972) in the Design of UCLA General 
Circulation Model. 



requirement that the fluxes of the surface layer and the mixed 
layer are consistent with each other. Similarly, surface 
humidity is determined by requiring equality of the surface 
layer moisture flux and the mixed>layer moisture flux. 

Since surface temperature is not known a priori , in 
the old parameterization, the surface temperature from the 
previous time step was us«d to calculate the new fluxes, which 
were then used to obtain t new surface temperature. The 
final surface temperature was then the average of the old and 
the new values. Finally, this surface temperature was used to 
recalculate the surface fluxes. In affect, it amounts to one 
cycle of iteration. Even with all this averaging, a 2-5t 
oscillation in the surface temperature occurred. This gave 
rise to corresponding oscillations in the surface fluxes. 

G OVERNING EQUATIONS 

In the GLAS model, a stable PEL forms whenever the 
ground is cooler than the air, and an unstable PEL forms whenever 
the ground is warmer than the air. The air temperature used in 
this comparison is the potential temperature at the lowest model 
level, level 9 of the GCM, which is nominally 945 mb. Of course, 
for saturated air a correction is necessary to account for the 
lapse rate modification by moisture. The parameterization 
proposed here is based on a unique solution for the surface 
temperature, which corrects the 2-6t oscillation found in the 
old parameterization. An "analytic" solution for Tg is obtained 
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for the stable PBL» and a curve-fit solution for Ts is obtained 
for the unstable PBL. The basic quantities and useful variables 
are defined below t 




ground temperature, R 

saturation mixing ratio at the ground 


Tg, Ogx physical and potential temperature of layer 9 
in the model 

qg • mixing ratio at level 9 of the atmosphere 

qs • surface mixing ratio 

U_ V- ■ surface U,V winds, in ms“^ 

9 # 9 


Ug^Vg ■ level 9 U,V winds, in ms 


-1 


w 


s 


“s* ♦ V,2 


K ■ Eddy transport coefficient of the mixed layer cal 
L “ latent heat of evaporation cal g"^ 

Cp B specific heat at constant pressure 
R » gas constant 
< » R/Cp 

e s ratio of molecular weights of water and air 

The temperature difference across the mixed layer, 58, is 
modified to reflect the moist adiabatic lapse rate as follows 
The lapse rate for dry atmosphere is obtained from 


dT 

37 


I'd 


For the saturated air, the relation is 


dT 

37 


^s= 




* 

h 9^ 

1 ♦ R 

L de 
1 + Cp dT 


m-lifl 


( 1 ) 


( 2 ) 
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If the relative humidity in the real atmosphere is rs, then rx>s 
is obtained by linear interpolation as follows: 

Fr “ (1 “ rs) F(j + rs Fg , (3) 

s 


The above calculation assumes that the saturated air mass in 
the PBL regime of a grid cell is proportional to the relative 
humidity of air involved. 



Equation (4) follows from (3) and (2) with the use of the following 
additional relations: 


* 

and eg = qgP/e 

Multiplying (4) by the boundary layer height, and using 


tfPg - Fj» AZp 

s 


(5a) 
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we obtain 


ATc • ATa + Yc • 

To solve the PBL equations to obtain the surface temperature, 
we define 


A 6 ■ Tg - 69 , 

«0 ■ Ts - 99 r 

and A0 - 60 ■ Tg - Ts » 

where 0g ■ Tg (pg/pg)* and 0_ s T_ and 0 gS T, 

(see Fig. 1). v v 


(7a) 

(7b) 

(7c) 


We solve for 60 for a given A0, then obtain Tg from (7c). 

The drag coefficient, (also equal to heat transport 
coefficient), is a linear function of surface geopotential 
over land and a linear function of surface wind over the oceans. 
Assuming that Cq is known for a given grid cell, the surface 
heat flux, Fg, may be obtained as follows: 


Fs “ P*Cp Dr ( 0g ~ 0g ) 


( 8 ) 


where 


Dr - CpWg3/{Wg2 - 7x(A0 - 60)} if Tg<Tg , (9a) 


or 


Dr 



(A0 - 60) } 


if Tg>Tg. (9b) 
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The mixed-layer heat flux, given by 


Pml - - pCp K/Zb(AT - ATc) , 

where K is the eddy diffusivity of the mixed layer. There are 
separate parameterixations for the calculation of K for the 
stable and unstable PBL's as follows: 

For the stable case 

K » 60. /(I + 40.*Ri) , 

where Bulk Richardson number r Ri» is defined as 

Ri = - «0g2b/l0g* {(Ug-Ug)^ + (V9 - Vg)2}]. 

For the unstable case 

K = 60. + 100. [1 - exp(-l,2 )1 • 

3Z 

By equating Fg and Fi-. , we find that 

DR(0g - 0s) = -K/Z 3 (AT - ATd " Yc ) “ K/Zb( «0 + Yc>- 
Here we have used AT - yTj - Yc “ (Tg - Tg) - (Tg - Tg)da " Yc 

» -l(Tg - Tg ) + Ycl 

da 


= -( 60 + Yc) • 

From (14), and with given values of 0g, Tg, 09 , Yc» Dr, and K, 
it is now possible to determine 50 and Og. The case for an 
unstable PBL is solved by a curve-fit relation between 


( 10 ) 


( 11 ) 


( 12 ) 


(13) 


(14) 
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66 ■ £(66, Cq, Hg). cado for a stable PBL is determined 

exactly by solving a quadratic equation in 66. 

(a) Solution for 66 (Stable Case) 

A stable case is obtained by using Equations (14), (9a), 
and (11). 


66 


{Ci,Ws^/lWs^ - 7 (66 - 66))}66 


CjjW VIW 32 - 7(66 - 66)1 +1 


(60/(1 - 40 


66gZ3/0gDw2)] 


where . 0^)2 + . Vg)2 


which rearranges to 

^(7 X eO/ZgCoWg^ - 4OgZ3/0gE^2j^j02 + 

— V ■ ^ 

A 

^1 + 6O/Z3C3W3 - 7 x 60 ./ZgCj)Wg^ + 40 gZ 3 Cp/ 6 gD„^)^ 66 « 66 
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from which 66 may be obtained. If A approaches zero in the 
quadratic equation A662 + b 66 ■ 66, the limiting solution is 

66 ■ - 66/B(l - A. 66/b2). (16) 

In the above formulation, the boundary layer height Zg is of 
the order of 500 m. However, a stable PBL is generally shallow. 
Its height is of the order of 100 m. In order to reflect the 
difference in height, the constant *40' is changed to '8.* 
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A limiting value of K ■ 2.0 is reached when the limiting 
critical Richardson Number 3.05 is attained (Deardorff, 1967). 
Accordingly, the minimum value of <6 is obtained combining 
(14), (9a) and (11) to give 

«0 - + (K„in/ZBCjjWa)(l-7.Ae/Wg2)] . (17) 

(b) Solution for 46 (Unstable Case) 

In an unstable case, (14), (9b) and (13) must be solved. 
Since (13) is transcendental, an exact solution is not possible. 
Hence, an attempt is made to obtain 66 as a function of 46, Wg, 
and Cq, for a range of values. Figures (2a), (2b), and (2c) 
show 66 as a function of 46 for values of Cq between .001 and 
.005 and surface wind magnitudes of 2-12 ms*"^. In all the 
graphs, the lines are the curve-fit solutions to match the 
points which are exact calculations. Prom these solutions, 
graphs are obtained for 66/46 versus Cp, for various values 
of the winds (Figure 2d) and 66/46 versus surface wind for 
various values of Co (Figure 2e). 

A simple functional form to obtain curve-fit relation 
66 and 46, Cq and Wg can be derived. Obviously, 

66 « f(46, Cd, Wg). 

However, since 66 is linear with A6, a suitable functional form 
will be 

66/46 • f(Co. Wg). 

nut at Wg * 0, 6 6/46 / 0. Thcroforo, assume the functional lot m 
to be 

60/46 - fl(CD,Wg) 4 f2(Co). 
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Some preliminary calculations indicate a general form to be 

«e/A6 - + f 2 (Co>, (21) 

which was approximated by 

«e/ae - + AjCp, (22) 

where Ai and A 2 are arbitrary constants. Using the method of 
least-squares, the constants A^ and A 2 were found to be 
0.1382 and 13.67, respectively. Therefore, the final form of 
the relation between 66 and A6 is: 

«6 « (0.1382 X Wg X + 13.67 Cp)Ae. (23) 

(c) Wind Field Modification by Surface Drag 

The surface drag force acting on a grid cell produces 
a change in the momentum of the air at level 9, which, in the 
GLAS parameterization, is the only layer directly affected. 

Also, for Aralcawa Grid B (Ara)cawa, 1972), the wind fields are 
defined at the secondary points, whereas the drag force is 
calculated at the primary points. In the old code, first the 
wind fields were interpolated to primary points, then momentum 
deficit was calculated for the interpolated wind fields. This 
momentum deficit was then reinterpolated to secondary points. 

This bac)c-and-forth interpolation does not allow a direct 
coupling between wind velocity and momentum deficit. 

The effect of this interpolation ma/ be minimal if 
every field is smoothly distributed. However, these wind 
fields are not snooth. Instead, large gradients are found. 


9 



particularly In the event of a growing 2-6x oscillation. The 
scheme of interpolation described above may result in spurious 
and sometimes systematic momentum transfers between grid cells, 
thus feeding these oscillations. Besides, and most importantly, 
any 2-<x pattern, if present, will be unaffected by friction. 

The new boundary layer calculation partly, if not completely, 
eliminates such a 2-«x oscillation. In the new procedure, 
first a factor. Dr, is calculated at each primary point, as 
before. It is then interpolated linearly to secondary points. 

The momentum deficit is now calculated by multiplying this 
factor by winds at level 9. Mathematically, the old method of 
obtaining AV was; 

*Vsi,j " {AVpi^j + AVpi+i^j + AVpi^j«l + AVpi+i^j-l} 

( 24 ) 


whore 


, i , j 


- Ji 

AP 




D 

Pi»j 


P . . 

Pi»3 


V 

Pi* j 


At 


It is replaced by the calculation given belo.^; 

t 

First define ®Ppi,j as follows 
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( 25 ) 


( 26 ) 


Now 

9 

Dr,1,j . 0.25. l%i, 5 




( 27 ) 


and 


AV 9 


Si, j 


^i,j . 


(2b) 
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The symbols and indices used in the above equations are as follows 
Indices irj represent longitude and latitude. Suffix 'P* stards 
for variable defined at primary and 'S' for secondary points 
of the grid on Scheme B. Pg is the mass of air in sigma 
level-9 in millibars, p is the density in 10*1 gm/cm^. At 
is time step in seconds. The drag factor » is defined by 
(9a) and (9b) for a stable and an unstable PBL, respectively. 

t 

Variable D^s is defined in Equation (28). Also, wind velocity 
updates are saved in an array. These are made simultaneously 
at the end. 

ANALYSIS OF RESULTS 

(1) First a time step invoking one call each of 
radiation and physics and three calls of hydrodynamics was 
completed to compare the results. The only striking difference 
was larger evaporative tendencies over oceans, and some large 
differences in surface temperature were noticed. Figure 3 shows 
the digitized maps of surface temperature differences. A 
large difference of 10-20*C may be seen in regions marked (I). 
These differences occur because in this region the air above 
the surface is very cold. The old model makes Tg close to 
the ground temperature, whereas the new model makes it close to 
the level-9 potential temperature. This is so because in the new 
calculation the eddy exchange coefficient has increased by an 
order of magnitude for the unstable atmosphere. The effect of 
this is to increase the forcing gradient of ~he surface fluxes. 

The larger temperature gradient, particularly over oceans. 
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increases Dj^ in equation (9b) and surface flux in equation (8). 
The same line of reasoning applies to moisture. However, over 
the land, increased surface fluxes will reduce the diurnal 
temperature oscillation. This suppresses evaporation in favor 
of increased sensible heat flux, because the diurnal oscillation 
of surface saturation humidity is several times the magnitude 
of the temperature oscillation. Thus, the net flux of moisture 
is reduced relative to the sensible heat flux. 

(2) In a 1-day simulation, the 2-6t oscillation was 
eliminated. Figures (4a) and (4b) show the changes in the evapo 
ration and sensible heat flux over land and ocean in the two 
hemispheres separately. The previous runs show a very noisy 
field compared to the new run. Consistent changes in daily 
averages of sensible heat flux and evaporation on individual 
grid points are now simulated. 

SUMMARY 


A new PBL parameterization has been tested with an 
exact solution for surface temperature instead of an iterative 
and time-averaged solution. Also incorporated is the original 
formulation of the eddy diffusion coefficient of Deardorff. 

The results from the new parameterization show desirable effects 
of increased evaporation over ocean and reduced evaporation over 
land. The proposed procedure of calculating surface temperature 
also eliminates the 2-fit oscillation in the surface fluxes of 
the old model. A short-range forecast revealed snail but 
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b«n«£lclal effects on surface temperatures t sea level pressursf 
and geopotential heights at 500 mb. 

ACKNOWLEDGMENTS 

We would like to thank Dr. J. Shukla for encouraging 
us to independently solve this problem. Dr. D. Randall suggested 
the technique for wind field modification by surface drag. We 
sincerely thank him for his input. Also, we gratefully acknowledge 
support and interest expressed by Dr. M. Halem. 


REFERENCES 


Arakawa, A., 1972t Design of the UCLA atmospheric general 

circulation model. Tech. Report No. 7, Dept, of Meteor., 

Univ. of Calif, at Los Angeles. 

Charney, J. G., W. J. Quirk, S. H. Chow, and J. Kornfield, 1977: 

A comparative study of the effects of albedo change on 
drought in semi-arid regions. J. Atmos. Sci. , 34 , 1366-138'^. 

Deardorff, J. W., 1967: Empirical dependence of the eddy 

coefficients for heat upon stability above the lowest 50 m. 

J. Appl. Meteor. , 6, 631-643. 

Halem, M., J. Shukla, Y. Mints, M. L. C. Wu, R. Godbole, G. Herman, 
Y. Sud, 1979: Comparisons of observed seasonal climate 

features with a winter and summer numerical simulation 
produced with the GLAS general circulation model. Report of 
the JOC Study Conf. on Climate Models: Performance, Inter- 

comparison and Sensitivity Studies, GARP Publ. Series No. 22 , 
207-253, WMO, Geneva, Switzerland. 

Somerville, R. C. J., P. H. Stone, M. Halem, J. C. Hansen, 

J. S. Hogan, L. M. Druyan, G. Russell, A. A. Lacis, W. J. 
Quirk, and J. Tenenbaum, 1974: The GISS model of the 
global atmosphere. J. Atmos. Sci. , 31, 84-117. 


13 




- parameterization 

PR SQ Ys^ 


o 

o 

o 

tl 


-• % 


J I 



•a c 
c 0 

•H -i-t 
3S4J -P 
U 3 

•M lOr-l 

0X0 

U (0 

n 

O 4J 
P ‘-H 

10 jX) 

>a< (I) 

> 

(OOP 

p-l p 

oja u 

•H 10 

P -P c 
to CO o 
> c 
P c 
u 3 

O P o 

*M OjC 
«P D) 

<D 

<<H m 
O 4J 

• o c 

CO . -H 

>o 0 

ir a 

CD Q 

*o o 0) 

p • 
Di'O (0 (0 

c c o 
•H (0 to c 

S C -H 

O to o 

43 <U 

CO'S 44 CO 

3 (0 to 

CO 44 •-( 
43-H 3 C 

ac u » 

^0 O^fH <0 
M <0 ^ 

o E O >0 




Figure 2(b). Same as Figure 2(a) except C-s0.003 for unstable PBL 



Figvire 2(c). Same as Figure 2(a) except Cj^ssQ.OOS for unstable PCO^ 
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Figure 2(d). Graphs showing 5d/d6 vs. Cq for various wind 

magnitudes (for unstable PBL) . Exact calculations 
are points shown on curve fit solutions drawn as 
lines. 



Figure 2(e). Graphs showing 66/AO vs. wind magnitude 
for various values of (unstable PBL) . 





